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I.  INTRODUCTION 


Comparing  two  high  resolution  models  of  combat  has  proven  to  be  a  difficult 
task.  This  thesis  uses  a  simple  discrete  dynamical  model  of  combat  as  a  surrogate 
for  the  complex  models  being  studied.  Our  surrogate  model  is  fit  to  simulated  battle 
data  generated  by  a  mathematical  model,  and  used  to  generate  probability  surfaces. 
These  surfaces  are  compared  using  a  randomization  test  [Ref.  1]  which  proves 
effective  at  identifying  battle  data  sets  from  similar  simulations  models  and  at 
distinguishing  between  battle  data  sets  from  different  simulation  models. 

Land  combat  has  evolved  into  a  series  of  complex,  combined  arms  battles  fought 
with  extremely  expensive  high-technology  weapons.  Many  recent  efforts  to  model 
this  process  have  relied  on  computer  simulations.  These  simulations  attempt  to 
capture  the  processes  of  combat  by  simulating  the  detailed  interactions  between 
individual  combat  elements  (e.g.  tanks,  armored  personnel  carriers,  and  artillery 
pieces).  As  a  result,  these  high  resolution  simulations  are  highly  complex  and 
expensive  in  their  own  right. 

Having  paid  the  price  to  develop  these  simulations,  users  rightly  want  to  know 
if  they  replicate  the  combat  processes  modeled.  This  question  has  to  a  large  extent 
gone  unanswered.  A  great  deal  of  data  analysis,  experimentation,  and  field  testing 
goes  into  developing  procedures  used  to  simulate  the  various  "microscopic"  combat 
processes,  such  as  searching  for  targets,  identifying  targets  as  hostile,  and  engaging 
hostile  targets.  Verifying  that  these  procedures  adequately  model  the  processes  for 
which  they  are  designed  is  relatively  straight  forward.  The  interactions  amongst  all 
the  various  underlying  processes  are  not  well  understood,  and  that  makes  the 
validation  of  the  overall  model  difficult. 

To  validate  a  high  resolution  combat  simulation  such  as  Janus(A)  [Ref.  2],  an 
army  battalion  task  force  air-land  combat  simulation,  one  might  try  comparing  some 
measure  of  effectiveness  (MOE)  from  an  actual  combat  engagement  to  the  same 
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measure  of  effectiveness  for  a  simulation  of  that  battle  in  Janus(A).  However,  not 
much  data  are  available  on  actual  engagements  of  high-tech  armored  forces,  and 
actual  combat  data  tends  to  be  clouded  by  haphazard  data  reporting  and  collection 
techniques.  The  next  best  procedure  might  be  comparison  to  data  collected  from 
realistic  field  training  exercises  such  as  those  conducted  at  the  United  States  Army’s 
National  Training  Center  (NTC)  at  Fort  Irwin,  California.  The  data  from  NTC, 
although  collected  more  systematically  than  combat  data,  also  has  its  shortcomings 
[Ref.  3],  and  it  is  expensive  to  collect.  For  these  reasons  it  is  difficult  to  get  the 
desired  number  of  replications  of  similar  battles  from  NTC. 

How  then  can  we  validate  a  model  such  as  Janus(A)  if  there  is  no  good  source 
of  data  for  comparison?  This  thesis  attacks  that  problem  by  exploring  ways  to  use 
analytical  (mathematical)  surrogate  models  of  high  resolution  simulations,  such  as 
Janus(A)  and  field  exercises  at  NTC,  fitted  to  small  battle  data  sets.  Chapter  II 
outlines  the  methodology  proposed  to  compare  two  combat  models.  Chapter  III 
examines  the  utility  of  simulated  annealing,  steepest  descent,  non-linear  optimization, 
and  multiple  regression  techniques  for  fitting  the  parameters  of  an  analytical 
surrogate  model  to  the  simulation  model’s  data.  Chapter  IV  presents  an  example  of 
the  methodology  to  compare  similar  and  dissimilar  battle  data  sets  generated  from 
a  mathematical  model,  and  finally  Chapter  V  concludes  this  study  with  a  discussion 
of  the  utility  of  the  methodology  and  proposals  for  further  research. 
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II.  METHODOLOGY 


A.  INTRODUCTION 

This  chapter  examines  the  methodology  developed  for  comparing  two  high 
resolution  models  of  combat  (referred  to  as  simulation  models  throughout  this 
thesis).  First,  we  examine  the  assumptions  necessary  for  our  method,  and  the 
surrogate  models  considered  (including  a  development  of  a  discrete  Lanchester 
model).  We  then  discuss  the  techniques  investigated  for  fitting  battle  data  to  the 
surrogate  model  and  the  measures  considered  for  comparison.  Finally,  we  examine 
the  actual  procedure  used  for  comparing  the  two  simulation  models  (i.e.  two  machine 
generated  battle  data  sets  representing  simulation  models). 

Throughout  this  chapter  we  will  refer  to  two  critical  items  which  are  defined 
below. 

•  Battle  trajectory.  A  battle  trajectory  is  a  numerical  record  of  battle  expressed 
as  the  number  of  each  weapon  system  still  active  in  the  battle  at  the  end  of 
each  time  increment. 

•  Battle  data  set.  A  battle  data  set  is  a  collection  of  battle  trajectories  from  the 
same  simulation  model. 

B.  ASSUMPTIONS 

Certain  assumptions  are  necessary  to  implement  this  methodology. 

1.  Battle  trajectory  data  are  obtained  as  the  raw  number  of  each  weapon 
system  actively  involved  in  the  battle,  recorded  at  the  end  of  each  time  interval,  as 
shown  in  Table  1;  all  time  intervals  are  of  the  equal  length. 

2.  It  is  assumed  that  each  weapon  system  involved  in  the  battle  is  capable  of 
engaging  any  other  weapon  system  involved  in  the  battle.  This  assumption  is 
required  when  using  Lanchester-type  equations  as  the  surrogate  model  (Lanchester 
Models  of  combat  are  discussed  in  Section  D). 
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TABLE  1.  SAMPLE  BATTLE  TRAJECTORY. 


Battle  Trajectory  from  NTC 


Time 

RTANK 

RAPC 

BTANK 

BAP( 

0 

60 

80 

50 

50 

5 

56 

73 

47 

46 

10 

52 

67 

45 

42 

15 

48 

61 

43 

39 

20 

45 

55 

41 

36 

25 

42 

50 

39 

33 

30 

39 

45 

37 

30 

35 

36 

41 

35 

28 

40 

34 

37 

34 

26 

45 

32 

33 

33 

24 

50 

30 

29 

32 

22 

55 

28 

26 

31 

20 

60 

26 

23 

30 

18 

BTOW 


3.  The  battle  trajectories  are  assumed  to  be  non-increasing  (i.e.  no 
reinforcements  are  allowed).  This  is  often  the  case  in  computer  simulation  models. 
This  assumption  forces  the  combat  processes’  mathematical  epresentation  to  be 
convex,  thus  simplifying  the  fitting  of  the  analytical  model. 

4.  Individual  trajectories  are  assumed  to  be  independent  of  all  other 

battle  trajectories.  This  is  not  a  problem  with  computer  simulations  that  are 
restarted  using  different  random  number  seeds.  However,  this  can  be  a  problem  with 
data  obtained  from  field  exercises  such  as  those  conducted  at  the  NTC.  Similar 
battles  conducted  by  the  same  unit  (or,  even  with  the  same  opposing  force)  can 
display  learning  curve  trends  which  could  result  in  correlation  between  different 
battle  trajectories. 

C.  THE  SURROGATE  MODEL 

Choosing  the  surrogate  model  is  a  critical  step  in  the  analysis  process.  The 
model  chosen  must  represent  the  underlying  processes  of  interest.  Ideally,  the  model 
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should  also  be  simple  enough  (i.e.  have  a  small  number  of  parameters  to  be 
estimated)  to  allow  good  fits  to  small  data  sets. 

The  model  chosen  for  our  initial  trials  was  a  discrete  analog  of  the  Lanchester 
aimed  fire  model.  This  model  was  chosen  because  it  is  an  attrition  model;  because 
it  has  a  small  number  of  parameters  to  be  estimated;  and  because  the  weapon 
systems  considered  were  all  direct  fire  weapons  (i.e.  tanks,  armored  personnel 
carriers,  and  anti-tank  weapons).  A  short  introduction  to  Lanchester-type  models  is 
provided  in  the  next  section. 

D.  IANCHESTER  MODELS  OF  COMBAT 

In  the  late  eighteenth  century  Antoine-Henri  Jomini  wrote  extensively  on  what 
he  called  the  principles  of  war  [Ref.  4].  One  of  these  principles  was  that  an  army 
should  mass  (concentrate)  its  forces  when  attacking.  Over  a  century  later,  in  1914, 
a  British  engineer  and  inventor  F.  W.  Lanchester  set  out  to  prove  a  hypothesis 
similar  to  Jomini’s  [Ref.  5].  Lanchester  conjectured  that  in  "modern  warfare"  the 
concentration  of  forces  was  an  effective  tactic.  To  prove  his  hypothesis  Lanchester 
developed  a  mathematical  model  of  "modern  warfare."  Lanchester  argued  that  with 
modern  weapons  it  was  possible  for  any  combat  element  to  engage  any  other  on  the 
battlefield,  thus  the  attrition  of  a  force  should  be  proportional  to  the  size  of  the 
opposing  force.  Lanchester  presented  the  following  model: 


dx 

dt 

dY 

dt 


-a  Y ; 

a  >  0  , 

*(0)  =  Xq  , 

(1) 

-bX; 

b>  0  , 

Y(0)  =  y0  , 

(2) 

where 

X(t)  =  X  force’s  strength  at  time  t, 

Y(t)  =  Y  force’s  strength  at  time  t  [Ref.  6]. 
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This  model  can  be  solved  by  dividing  equation  (1)  by  equation 

(2): 


dx 

dt 

-a  Y 

(3) 

dY 

-bX  ' 

dt 

dX 

-a  Y 

(4) 

dY 

-bX  ' 

XdX  = 

-a  YdY. 

(5) 

We  then  integrate  both  sides: 


t  t 

J  -bxdx  =  [  ~a  YdY, 

0  0 

Z  0  Z  o 

-b{X2{t)  -  X2(0))  =  -a(Y2(t)  -  Y2(  0)).  (8) 

Thus  we  have: 

b(x02  -  X2  ( t)  )  =  a( y02  -  Y2  ( t) ) ,  (9) 

which  is  known  as  Lanchester’s  Square  Law. 

Lanchester’s  model  is  a  continuous  approximation  of  the  discrete  combat  attrition 
process.  This  formulation  presents  two  problems  for  direct  use  as  a  surrogate  model. 
First,  in  order  to  use  Lanchester’s  model  to  simulate  combat  on  a  computer,  it  is 
necessary  to  approximate  the  differential  equations  with  difference  equations. 
Second,  the  observations  we  are  using  to  estimate  the  parameters  in  our  model  are 


(6) 

(7) 
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taken  at  discrete  time  intervals.  It  thus  makes  more  intuitive  sense  to  formulate  the 
model  with  difference  equations  from  the  start.  Lanchester’s  model  expressed  as  a 
system  of  difference  equations  is: 

■^n+l  Xn  =  “ a  Yn  !  <3  >  0  /  Y0  =  y0  ,  ( 10 ) 

YB.1  -  Yn  =  -bxn;  b>0,  X0  =  x0.  (11) 


By  iterating  equation  (10)  one  step  we  have: 

Xn* 2  ~  ^n*  1  =  (12) 

We  then  subtract  equation  (10)  from  equation  (12)  which  yields: 

*n+2  -  ^  Xn+1  *  Xn  =  -a(Yn+1  -  Yn),  (13) 

and  by  substitution  we  get, 

^2  -  2*ml  +  =  -a  (  ~bXn)  .  (14) 

Thus, 

Xn.2  =  2Xb.1  +  (ab  -  l)  x„.  (15) 


which  is  a  second  order  discrete  dynamical  system.  It  can  be  solved  using  standard 
techniques  for  second  order  difference  equations  [Ref.  7].  Its  characteristic  function 
is  given  by:  r2  -  2r  +  (1  -  ab)  =  0. 

This  equation  has  two  real  roots, 


Thus  the  general  form  of  the  solution  to  the  discrete  dynamical  system  for  the  X 
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force  is  given  by: 

Xn  =  C0(  1  +  /aB)n  +  q(l  -  /aj b)n ,  (17) 

where  C0  and  Q  are  constants  determined  by  the  initial  conditions.  Applying  our 
initial  conditions: 

Xo  =  Xo  and  Xj  -  Xo  =  -a  y0, 
we  obtain  the  following  two  equations  in  two  unknowns: 

X0  =  x0  =  C0(l  +  +  Cx(  1  -  yfaBf ,  (18) 

Xx  =  x0  -  a  y0  =  C0(  1  +  \J~aB)1  +  Cx(l  -  JaB)1  .  (19) 

Solving  for  C0  and  Cx  we  get: 


Thus  the  particular  form  of  the  solution  for  the  X  force  is: 

*„ =  -r  "  \  i  y»  (1  +  '^E)”  *  t  *  \  i  y°  (1  ■  ■  <20) 

and  similarly  for  the  Y  force: 

r„-  f  x,  (1  *  fad)"  *  -  \  |  X0  (l  -  JaBV  .  (21) 

The  graph  of  these  solutions  with  two  different  values  of  v^ab  is  shown  in  Figure  1. 
From  this  graph  we  see  that  the  solution  trajectories  are  approximately  linear  over 
the  domain  of  interest  (the  first  quadrant)  with  slopes  determined  by  the  value  of 
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Discrete  Lancf^ter  wxxi 


6  B  ID  12 


Figure  1.  Solution  Trajectories. 


9 


Jab.  The  larger  the  value  of  Jab  the  more  intense  the  attrition  process  and  thus  the 
steeper  the  slope.  We  can  interpret  the  Jab  as  a  measure  of  the  pace  or  intensity 
of  battle  (see  Section  F). 

Another  value  of  interest  is  the  time  at  which  the  losing  force  would  be 
annihilated.  We  are  unable  to  calculate  this  value  for  the  continuous  Lanchester 
model  because  the  solution  trajectories  go  to  zero  asymptotically. 

For  the  discrete  model  the  annihilation  time  of  the  losing  force  is  calculated  by 
setting  X(nx)  =  0  and  solving  for  nx.  Thus  we  have: 


X(nx)  =  0  = 


2  ~  \ 


i  » 


(l  +  fa£)nx  + 


+  N 


~B  y° 


(l  -  yfa5)n* , 


N 


±  yn  - 

b  y°  2 


(l  +fab)nx  - 


+  N 


3yo 


(l  -  Jabfx , 


1  -  y lab 
[  1  +  sfah 


N 


-  -  — 
b  y°  2 


*  N 


iy° 


Now  provided  the  right  hand  side  of  this  equation  is  positive,  which  we  show  below, 
we  can  take  logarithms  on  both  sides: 


In 


nx  = 


M 


iy° 


*0  A 

a 

T  \ 

T  y° 

In 


•fab 


l  +  fab 


N 


i* 


-  ^S.  >  0  ,  Jah  <  1 
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Similarly  from  equation  (21)  we  get  the  annihilation  time  for  the  Y  force: 

We  no\.  show  that  the  quantities  in  the  above  equations  are  positive.  We  know  that 
if  \/ab  >  1  then  at  least  one  of  the  attrition  coefficients  a,b  >  1.  This  implies  that  in 
the  battle,  at  least  one  of  the  forces  would  be  annihilated  in  one  time  period.  Thus 
we  are  only  interested  in  the  case  when  ^ab  <  1,  and  we  have  adjusted  the  length 
of  our  time  periods  to  insure  that  a,b  <  1.  Thus  in  order  to  show  that  we  can  always 
determine  the  battle  termination  time  nx  we  need  only  show  the  following. 

Lemma. 


Assume  *0,  y0  >  0  and  0  <  a,b  <  1. 


Proof: 


Q.E.D. 
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Using  the  formulas  we  have  derived  for  the  battle  termination  time  we  can  now 
develop  conditions  for  determining  which  force  will  be  victorious.  Setting^  =  ny 
we  get: 


Multiplying  out  the  terms 
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that  is, 

y?_  =  b 
V  ~  ’ 

y  2  . 

Therefore  if  —  =  —  then  the  battle  is  a  stalemate, 

x?  a 
•*o 

y  2  t 

else  i if  —  >  —  then  the  Y  force  wins , 

r  2  a 

■*o 

y  2  . 

if  —  <  —  then  the  X  force  wins. 

V  2  ti 


This  analysis  of  our  model  suggests  a  number  of  natural  measures  of  battle  which 

we  will  consider  in  Section  F. 

1.  Non-Homogeneous  Lanchester 

Our  model  up  to  now  has  considered  only  homogeneous  forces,  such  as  a 
battle  of  tanks  versus  tanks.  The  models  we  wish  to  compare  contain  not  only  tanks, 
but  armored  personnel  carriers  and  anti-tank  weapons  (they  may  also  contain  indirect 
fire  weapons  such  a  artillery;  however,  we  will  not  consider  them  here).  To  include 
the  differences  in  the  attrition  potential  of  the  various  weapon  systems,  we  assign 
each  weapon  a  weight,  which  we  call  its  "fire  power  index."  The  total  attrition 
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potential  of  the  X  and  Y  forces  is  calculated  as  a  linear  combination  of  the  individual 
fire  power  indices  [Ref.  8j. 


I 

*  *  Ll  *j  ■  cm 

j 


Where; 

Y  *  the  total  fire  power  of  the  Y  force, 
y,  =  the  number  of  Y  systems  of  type  i, 
and 

fy  s  the  fire  power  index  of  a  Y  system  of  type  i, 

and  similarity  for  the  X  force. 

2.  Stochastic  Lanchester 

Our  model  has  been  deterministic  thus  far,  with  the  battle  results  being  fixed 
by  the  initial  force  levels  on  each  side  and  the  fixed  attrition  coefficients.  There  is 
intuitively  much  more  to  who  wins  a  battle  than  simply  the  initial  force  levels. 
Leadership,  training,  momentum,  and  terrain  are  just  a  few  of  the  critical  factors  in 
determining  the  outcome  of  any  real  battle.  Introduction  of  these  factors  explicitly 
into  our  model,  even  if  feasible,  would  make  it  too  complex  to  fit  to  small  data  sets. 
We  take  a  simpler  approach  by  adding  these  factors  as  stochastic  noise.  Gang  back 
to  the  difference  equation  models,  let  us  examine  how  this  could  be  done.  One 
simple  approach  is  to  add  a  random  term  to  each  equation: 

-  X,  -  -aYn  +  Zx(n) , 

y*.  1  -  Yn  =  -bXn+Zy{n),  (23) 

where  Zx  and  Zy  are  random  variables  with  given  distributions  whose  parameters  are 
to  be  estimated. 
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This  approach,  however,  suggests  that  the  noise  is  independent  of  the  force 
sizes  involved.  This  is  intuitively  unappealing  as  history  shows  the  more  complex  the 
battle  the  more  confusion  likely  to  be  present.  We  can  incorporate  this  idea  into  our 
model  by  postulating  that  some  factors,  such  as  combat  readiness  and  momentum, 
may  be  dependent  on  force  size;  while  other  factors,  such  as  leadership  and  training 
may  be  independent  of  force  size.  Thus  we  separate  these  factors  into  two  random 
terms: 

=  ~A(n)Yn+Zx(n), 

-  Yn  =  -  B(n)  Xn  +  Zy(n) ,  (24) 

where  A(n)  and  B(n)  are  random  variables  representing  the  force  size  dependent 
factors,  and  Zx(n)  and  Zy(n)  are  random  variables  representing  the  force  size 
independent  factors. 

Further  we  assume  that  these  factors  are  predetermined  at  some  level  at  the 
beginning  of  each  battle  based  on  the  training  and  maintenance  preparation  of  the 
forces  involved.  These  values  then  vary  about  these  fixed  values  based  on  noise 
introduced  by  the  confusion  and  stress  of  combat.  We  therefore  postulate  that  for 
any  given  battle  the  random  variables  (A(n),  Zx(n),  B(n),  Zy(n))  can  be  represented 
by  a  constant  plus  some  error  (noise)  random  variable  which  is  independent, 
identically  distributed  normal  with  mean  zero  and  variance  o2. 

Thus  we  write: 

-A{n)  =  -a  +  eA(n), 
z>)  =  z*  +  ez,(«). 

-Bin)  =  -  b  +  eBin), 

2yin)  =  zy  +  ezin).  (25) 
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Substituting  these  expressions  into  equation  24  yields  the  following  system  of 
equations: 

"  Xn  =  [-<*  +  e»]  Yn  +  [Zx  +  ez/n)l  » 

=  [-*  +  ^n)]Xn  *  [zy  +  ez(n)].  (26) 

It  should  be  noted  that  in  these  equations  the  a,  zx,  b,  and  Zy  are  only  constant  over 
the  duration  of  battle.  They  will  vary  from  battle  to  battle  based  on  the  combat 
preparation,  force  mix  (e.g.  the  ratio  of  tanks  to  APCs),  and  the  commander’s 
concept  of  operations.  Thus  over  the  battle  data  set  we  will  have  a  set  of  values 
{ai»  zv  bp  zy,}  f°r  each  battle  which  can  be  considered  to  be  realizations  of  the 

random  variables  A,  Zx,  B,  and  Zy,  where  A,  Zx,  B,  and  Zy  are  jointly  distributed  with 
parameters  to  be  estimated. 

Finally,  we  include  the  individual  attrition  potential  of  each  weapon  system 
(fire  power  indices)  from  equation  22. 

This  completely  specifies  our  model  as  follows: 

-  [-«  +  «»]  Yn  +  K  +  Gz(n)], 

Yn. i  -  F  =  [  b  +  eB(n)]Xn  +  [zy  *  ez(n)l  (27) 

where, 

X.  -  Y>-T.  f„h' 

i  i 

fx  =  the  fire  power  index  of  X  force  weapon  system  type  j, 

fy  =  the  fire  power  index  of  Y  force  weapon  system  type  j, 

Xj.  =  the  number  of  X  weapon  systems  of  type  j  active  in  the 
battle  at  the  end  of  time  increment  i, 
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Yjt  =  the  number  of  Y  weapon  systems  of  type  j  active  in  the 
battle  at  the  end  of  time  increment  i, 

eA(ri)  is  iid  NormaKO.o^2), 
ez  (n)  is  iid  Normal(0,oz  2), 
eB(n)  is  iid  Normal(0,o/), 
ez(n)  is  iid  Normal(0,o72), 

(note  eA(n),  ezjji),  eB(n),  and  ez(n )  are  also  independent  of  each  other) 

[a,  zx,  b,  zy]  is  a  realization  of  \A,  Zx,  B,  ZJ 
and  a,b  >  0. 

E.  FITTING  TECHNIQUES 

A  number  of  techniques  are  evaluated  in  this  thesis  for  estimating  the  parameters 
of  a  system  of  differential/difference  equations  from  small  battle  trajectories.  Using 
the  non-linear  programming  solver  MINOS  5.2  with  the  General  Algebraic  Modeling 
System  (GAMS)  to  do  a  simultaneous  least  squares  fit  on  the  data  was  found  to  be 
the  most  effective.  However,  suppose  one  is  willing  to  make  the  simplifying 
assumption  that  in  an  individual  battle  the  coefficient  random  variables  in  one 
equation  are  independent  of  the  coefficient  random  variable  in  the  other  equation 
(i.e  that  [A,  ZJ  and  [B,  ZJ  are  independent,  but  A  and  Zx  are  still  considered  to  be 
dependent,  as  are  B  and  Zy).  In  th;s  case  the  parameters  can  be  estimated  by 
standard  regression  techniques.  The  validity  of  this  assumption  can  be  tested  by 
fitting  a  small  number  of  battle  trajectories,  using  the  GAMS  model  which  considers 
each  equation  simultaneously,  and  comparing  the  resulting  estimated  coefficients  w ith 
estimates  obtained  from  the  standard  h'near  regression  fits  of  the  same  battle 
trajectories  (i.e.  those  obtained  by  considering  each  equation  separately).  If  the 
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comparison  is  favorable  (i.e.  the  difference  between  the  parameters  fit  by  the  two 
models  is  sma!1),  the  savings  in  the  computation  f'me  required  for  parameter 
estimation  is  great. 

F.  MEASURES  OF  BATTLE 

In  this  section  we  consider  measures  derived  from  the  difference  equation  model 
with  only  force  size  dependent  terms;  that  is,  all  factors  are  considered  to  be 
dependent  on  force  size.  With  this  assumption  equation  27  becomes: 

^(n)]  Yn  , 

Y..i  ~Yn  =  [~b  +  eB(n)]  Xn.  (28) 


In  defining  the  measures  of  battle  the  following  definitions  will  be  used. 

AXn  =  the  change  in  the  fire  power  of  the  X  force  during 
the  time  increment  [n,  n+ 1];  so  AXn  =  Xn+1  -  Xn. 


AYn  s  the  change  in  the  fire  power  of  the  Y  force  during 
the  time  increment  [n,  n+1];  so  AYn  =  Yn  +  1  -  Yn. 

Xn  s  fire  power  of  the  X  force  at  the  beginning  of  the 
time  increment  [n,  n+  1]. 

Yn  =  fire  power  of  the  Y  force  at  the  beginning  of  the 
time  increment  [n,  n+  1]. 

N  s  the  number  of  time  increments  in  the  battle 
trajectory. 


a  =  estimate  of  a  for  a  single  battle; 


a  = 


i=0 


N 


y-i 


A  Y 


b  =  estimate  of  b  for  a  single  battle;  ,  X 

b  = 


N 
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1.  Pace  of  Battle 

The  pace  of  battle  is  defined  to  be:  \jab  .  When  our  discrete  model  is  fit  to 

an  entire  battle  trajectory  the  pace  of  battle  provides  an  overall  measure  of  the 
violence  with  which  the  battle  was  fought.  It  may  be  of  even  greater  value  when 
calculated  at  each  individual  time  increment  [n,n  + 1]  of  the  battle  thereby  providing 
a  picture  of  how  the  intensity  of  battle  varied  over  time. 

2.  Relative  Combat  Power 

The  relative  fire  power  is  defined  to  be:  When  our  model  is  fit  to  an 

a 


entire  battle  trajectory,  calculation  of  the  relative  combat  power  gives  the  average 
relative  combat  power  of  the  two  opposing  forces.  Calculation  of  the  relative  combat 
power  during  each  time  increment  [n,  n+  1]  again  provides  a  picture  of  how  relative 
combat  power  varies  over  time.  This  can  be  viewed  as  how  the  tide  of  battle 

fluctuates  with  time. 

3.  Probability  of  Victory 

Based  on  our  model  (equations  28),  the  probability  of  victory  for  the  X  force 


can  be  approximated  by:  P[X  Wins ]  -  P  — <  —  .  Note  that  this  is  only  an 

y  2  A 
o 


approximation  since  we  are  ignoring  the  noise  terms  in  equation  28.  This  measure 
allows  us  to  measure  how  the  distribution  of  the  relative  combat  power  and  initial 
force  ratios  impact  the  battle  outcome. 
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4.  Elasticity 

Our  final  measure,  although  not  coming  directly  from  our  model,  is  composed 
of  the  same  data  elements  as  the  previous  measures.  It  is  known  as  the  elasticity 
[Ref.  9]: 

*Xn  K 

^  =  =  Xn  ' 

Yn 


Calculating  elasticity  at  each  time  interval,  we  obtain  a  picture  of  which  force  is 
winning  or  which  force  has  the  momentum.  Elasticity  is  similar  to  relative  combat 
power  in  that  it  also  provides  a  picture  of  how  the  tide  of  battle  varies  over  time. 

5.  Best  Measures 

Of  these  four  measures,  using  a  temporal  trace  of  the  pace  of  battle  in 
conjunction  with  the  relative  combat  power  provides  the  most  descriptive  graphical 
representation  of  each  battle  process;  providing  a  view  of  how  the  violence  and  tide 
of  battle  plays  out  over  time.  The  probability  of  victory  appears  to  provide  us  with 
a  useful  measure  for  analytical  comparison. 

Considering  possible  ways  of  computing  the  pace  of  battle  and  the  relative 
combat  power  suggests  two  methods  of  comparison. 

1.  If  we  calculate  each  measure  at  each  time  interval  then  we  get  a 
temporal  trace  of  the  pace  and  relative  combat  power  during  the  battle  (see  Figure 
2).  To  compare  these  traces  for  two  different  battle  trajectories  we  measure  the 
distance  between  them  (see  Barr,  Weir,  Hoffman  [Ref.  10]). 
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Figure  2  Temporal  trace  of  pace  and  relative  combat  power. 
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2. 


If  we  calculate  the  average  relative  combat  power 


A 

E 


B 

A 


N- 1 


E 


N 
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where  i  =  index  number  of  the  battle  trajectory, 

N  =  number  of  battle  trajectories  in  the  battle  data  set, 

over  each  battle  we  might  be  able  to  compare  these  measures  directly  using  the 
standard  t-test  since  the  distribution  of  relative  combat  power  appears  empirically  to 
be  normally  distributed. 

A  methodology  for  using  the  probability  of  victory  is  described  in  detail  in  the 
next  section. 

G.  COMPARING  THE  HIGH  RESOLUTION  MODELS 
1.  Methods  of  Comparison 

Regardless  of  the  fitting  technique  used  for  estimating  the  parameters  of  the 
model  (see  Chapter  III),  the  process  of  fitting  the  proposed  analytical  models  to  the 
simulation  data  results  in  a  set  of  estimates  of  the  realizations  of  the  coefficient 
random  vectors.  We  desire  to  use  the  discrete  model  with  the  estimated  coefficients 
to  compare  the  high  resolution  simulation  models.  To  do  this,  we  need  to  identify 
measures  that  can  be  used  to  detect  differences  or  similarities  in  the  underlying 
structure  of  the  high  resolution  models.  Direct  comparison  of  the  coefficient  random 
vectors  can  rely  on  independence  between  the  coefficients  and  on  sufficiently  large 
data  sets.  To  avoid  having  to  make  an  undesirable  independence  assumption,  and 
to  utilize  smaller,  less  costly  data  sets,  we  examined  two  indirect  approaches  for 
comparing  our  simulation  models. 

(1)  Fighting  simulated  battles.  This  involves  choosing  coefficients  from 
the  empirical  distributions  of  the  fitted  coefficients,  and  using  them  to  run  numerous 
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replications  of  a  simulated  battle.  The  MOEs  from  these  simulated  battles  are  then 
compared. 

(2)  Comparing  Probability  Surfaces.  Using  the  estimated  coefficient’s 
empirical  distribution  we  generate  probability-of-win  surfaces  for  various  initial  force 
ratios.  Surfaces  generated  from  two  different  battle  data  sets  are  compared  to  each 
other  numerically  to  obtain  a  difference  measure.  A  randomization  test  is  then  used 
to  compare  the  two  battle  data  sets. 

2.  Fighting  Simulated  Battles 

Our  first  method  uses  the  fitted  analytical  models  to  generate  larger  data 
samples  for  comparing  the  MOEs  from  each  model  being  considered. 

We  draw  pseudo  coefficients  from  empirical  distributions  of  the  estimated 
coefficient  random  variable  realizations.  These  coefficients  are  then  substituted  into 
the  surrogate  analytical  models  and  a  single  battle  is  fought  using  this  model,  yielding 
a  battle  trajectory.  We  then  draw  a  second  set  of  coefficients  and  repeat  the  process. 
This  continues  for  a  predetermined  number  of  replications. 

From  each  simulated  battle  trajectory  generated  in  this  fashion  we  observe 
the  MOE  of  interest.  The  empirical  distributions  of  the  MOEs  corresponding  to  the 
models  of  the  two  simulation  models  are  compared  directly. 

One  difficulty  with  this  method  is  determining  the  number  of  replications 
required  to  distinguish  between  two  different  battle  data  sets  from  two  different 
simulation  models  while  still  being  able  to  recognize  when  two  battle  data  sets  are 
from  the  same  simulation  model.  We  attack  this  problem  by  computing  the  sample 
means  (p,  and  p2),  and  sample  variances  (d,2  and  o22)  of  the  MOE  estimated 

directly  with  the  data  from  the  simulation  model.  This  information  is  then  used  to 
estimate  the  number  of  observations  required  to  obtain  a  confidence  interval  of 
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with  degrees  of  freedom  associated  with  . 

This  method  of  comparison  tended  to  transform  (or  shift)  the  mean  value  of 
the  MOE  upward,  and  to  introduce  additional  noise.  That  is,  the  variance  of  the 
MOE  from  the  simulated  (fitted)  system  was  greater  than  the  variance  of  the  MOE 
in  the  original  sample.  Since  our  intent  was  to  tighten  the  confidence  interval  on  our 
estimated  MOE,  the  introduction  of  additional  noise  was  counterproductive,  making 
the  comparison  less  sensitive.  Thus,  although  somewhat  effective  for  identifying 
similarities  and  differences  in  the  underlying  battle  structure,  it  was  determined  that 
this  approach  should  not  be  pursued  further. 

3.  Comparing  Probability  Surfaces 

From  our  analysis  in  Section  D,  we  know  that  if  (Yo)2/^)2  <  b/a  then  the 
X  force  will  be  victorious.  Since  we  are  assuming  that  A  and  B  are  random 
variables,  then  R  =  B/A  will  also  be  a  random  variable.  We  can  use  observations 
of  R  obtained  from  our  fitting  process  to  build  an  empirical  cumulative  distribution 
function  (cdf)  p  for  R.  This  empirical  cdf  determines  a  probability  surface  FR  over 
the  domain  {(xo,y0)  :  x,<Xo<xu  ,  y,<y0<yu  } 

x  2 

where:  FR{x0,y0)  =  P  R  <— 

>’o 

Xq  s  the  initial  fire  power  of  the  X  force, 
x,  =  the  lower  bound  on  the  initial  X  fire  power  values. 
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xu  s  the  upper  bound  on  the  initial  X  fire  power  values, 
and  similarly  for  the  Y  values. 

Suppose  we  wish  to  compare  two  sets  V  and  W  of  battle  data  (e.g.  one 
consisting  of  several  battle  trajectories  from  Janus(A)  and  one  consisting  of  several 
NTC  battle  trajectories).  A  natural  measure  for  this  comparison  is  the  "volume" 
between  the  two  surfaces: 


y  u 

v(W  =  £  E  |  ffl/Vo)  -  For(Wo>  I  A* Ay, 

>«  *  y/  *o m  *i 

where  Ax,  Ay  define  the  resolution  of  the  partition  of  the  domain  of  FR  (Ax  =  Ay  = 
1  was  used  in  our  trials). 

A  question  is,  "What  does  a  given  value  of  u(V,W)  mean?"  For  example, 
does  a  value  of  800  indicate  that  the  battle  sets  are  the  same  or  different?  To 
address  this  question  we  first  compare  probability  surfaces  with  a  randomization  test 
as  follows. 

Given  our  surrogate  model  (equation  28)  and  two  battle  data  sets  V  and  W 
consisting  of  n  battle  trajectories  each:  V  =  {v,,  v2, ... ,  vn}  and  W  =  {w,,  w2, ... ,  wn}, 
where  v,  a  the  i,h  battle  trajectory  from  simulation  model  V,  and  w,  a  the  i,h  battle 
trajectory  from  simulation  model  W,  we  estimate  the  coefficients  a  and  b  for  each 

battle  trajectory  in  V  and  W.  We  denote  these  estimates  bya  ,  a  ,  b  ,  and  b  , 

vi  wi  vi  wt 


where  av  a  the  estimate  of  a  obtained  by  fitting  battle  trajectory  vjt  and  the  other 

estimates  are  defined  in  the  same  manner. 

We  define  the  sets  R^,  and  as  follows: 


R v  =  {V  v  •  •  • .  where  \  = 


and 


R 


W 


We  now  define  an  empirical  cumulative  distribution  function  fr  for  a  set  of 
observed  values  of  the  random  variable  R;  jr(1),  r(2),  •  •  *  ,  r(lI)J  where 

r(l)  *  r(2)  s  '  '  ‘  sr(»)  • 


max  |i  j  r(<)  $  r} 

n 


[Ref.  13]. 


Examining  the  Xo-yo  plane  determined  by  a  finite  range  of  values  for  the  initial  fire 
power  of  the  X  force  (X^  and  Y  force  (Y0),  that  is,  the  plane 

{ (Vo)  I  XI  *  x0  *  xu>  yi  *  yo  *yu }  where  xi>  x«’  yi.  yu  are  integer  lower  and  upper  bounds 

on  the  range  of  values  for  Xq  and  y0,  we  partitioned  this  plane  into  a  grid  of  square 
intervals 

|  (x,y)  |  xt  +  k  s  x  <  +  yt+k  <.  y  <  yt  +  (k+l)  J  ,  k  =  0,  1, ,  m-1;  where  m  = 


xu  -  xi  =  yu  -  yi- 

We  now  use  the  definition  of  j?  and  the  partition  on  the  domain  of  F  to 

K  K 

define  a  probability-of-win  surface  S.  Given  a  set  of  observations  of  the  random 


variable  R,  {rp  r2, 


'  ,  rN j,  of  the  ratios 


choose  a  random  subset  R’  of  R. 


We  then  define 


V  =  y0>  /)  I  xi  *  xo  <  xu>  yt  *  y0<  /  =  fr- 


ho2'' 


r  2 

V*o  j 


We  define  our  test  statistic  to  be: 

r. 


“(V5*;)*  £  E  I'Vr 


y<r>/  *o=  xi 


2  'j 

yo 

_ 

2 

2 

(xo) 

26 


V(-V  SK-)  is  the  volume  between  the  probability-of-win  surfaces  generated  by  the 
sets  /?*  and  J?2 ,  where  R’  and  /?2*  are  sets  of  observed  values  of  the  random  variable 
R. 

To  determine  if  the  two  simulation  models  V  and  W  are  different,  we  first 
compute  T  =  SRJj.  This  is  the  observed  value  of  our  test  statistic.  We  then 

form  R  =  Rv  U  Rw.  Let  P  be  the  set  of  all  possible  partitions  R  into  two  subsets  Rj 

and  R2  such  that  both  R,  and  R,  contain  n  elements  and  /?,  f)  =  0.  A  specific 

partition  p  e  P  can  be  built  by  selecting  n  elements  of  R  at  random  (without 
replacement)  as  the  set  R(  and  letting  the  set  R2  be  R\R,.  We  can  then  compute 

vjSv  S^J.  Our  randomization  distribution  is  defined  to  be  the  set  of  values 

v  for  all  p  e  P.  We  compare  our  observed  value  T  with  a  sample  from  the 

randomization  distribution,  obtained  by  randomly  choosing  M  partitions  p  e  P  and 
evaluating  v^,  S for  each  p.  If  the  value  of  T  is  in  the  right  tail  of  this  sampled 

distribution  (i.e.  an  extreme  value)  we  infer  that  V  and  W  are  from  different 
simulation  models.  And  if  T  is  in  the  left  tail  of  the  randomization  distribution  we 
infer  that  V  and  W  are  from  the  same  simulation  model. 

Since  the  number  of  values  in  the  randomization  distribution  may  be  as  high 

as  -ip”)  =  @n)[  jt  woui(]  be  impractical  to  compute  each  value.  However,  it  is 
2  V n  )  4 (n!) 

reasonable  to  take  a  sample  from  this  distribution  for  testing  the  observed  value  of 
T  [Ref.  1]. 

In  our  studies,  this  technique  proved  to  be  effective  (results  are  given  in 
Chapter  IV)  at  identifying  difference  in  battle  data  sets  V  and  W,  when  V  and  W 
came  from  different  generating  systems  (see  Chapter  IV ,  Section  A  for  a  description 
of  the  generating  systems  used),  while  still  recognizing  similarity,  when  V  and  W 
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were  generated  by  the  same  system.  This  technique  is  demonstrated  in  Chapter  IV. 
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III.  FITTING  TECHNIQUES 


A.  INTRODUCTION 

Our  approach  to  finding  an  appropriate  technique  for  fitting  the  coefficients  of 
a  system  of  simultaneous  difference  equations  has  been  experimental.  We  first 
examined  the  simulated  annealing  technique  suggested  by  Ingber  [Ref.  2].  We  then 
considered  methods  that  required  less  computation  time.  This  chapter  is  a 
chronological  account  of  our  research,  we  thus  present  our  results  in  that  order. 

B.  SIMULATED  ANNEALING 

The  simulated  annealing  algorithm  was  developed  by  Metropolis  in  1953  [Ref. 
1 1]  to  simulate  the  physical  annealing  process  studied  in  statistical  mechanics.  It  was 
first  suggested  as  a  technique  for  solving  combinatoric  optimization  problems  in  1983 
by  Kirkpatrick,  Gelatt,  and  Vecchi  [Ref.  12].  The  simulated  annealing  algorithm 
combines  local  optimization  techniques  with  a  random  walk  process  to  reduce  the 
chance  of  becoming  trapped  in  a  local  optimum.  The  algorithm  shown  in  Figure  3 
begins  with  an  initial  solution  and  generates  a  proposed  neighboring  solution  at 
random.  If  the  proposed  solution  is  better  than  the  current  solution  (i.e.  a  downhill 
move)  then  the  proposed  solution  is  accepted  as  the  new  current  solution  with 
probability  one.  If  the  proposed  solution  is  worse  than  the  current  solution  (i.e.  a 
uphill  move)  then  it  is  accepted  with  a  probability  based  on  the  magnitude  of  the 
uphill  move  and  the  current  value  of  the  control  parameter  (referred  to  as  the 
temperature  of  the  process).  Thus,  small  uphill  moves  are  more  likely  to  be  accepted 
than  large  ones,  and  all  uphill  moves  are  more  likely  to  be  accepted  at  higher 
temperatures  than  at  lower  temperatures.  The  process  is  run  for  a  large  number  of 
repetitions  at  each  temperature,  and  reductions  in  temperature  are  made  according 
to  some  "cooling  schedule."  Changes  to  the  initial  solution,  initial  temperature,  and 
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cooling  schedule  can  have  dramatic  effects  on  the  convergence  properties  of  the 
algorithm  and  most  often  must  be  arrived  at  by  experimentation. 

Simulated  annealing’s  main  utility  is  in  solving  hard  combinatorial  optimization 
problems  for  which: 

•  no  good  problem  specific  heuristic  algorithms  have  been  developed  (simulated 
annealing  does  not  compete  well  with  problem  specific  techniques,  like  those 
developed  for  the  traveling  salesman  problem,  see  Johnson,  Aragon,  McGeoch, 
and  Schevon  (1989)  [Ref.  14]); 


1.  Get  an  initial  solution  S. 

2.  Get  an  initial  temperature  T>0. 

3.  While  not  yet  frozen  do  the  following. 

3.1  Perform  the  following  loop  L  time. 

3.1.1  Pick  a  random  neighbor  S’  of  S. 

3.1.2  Let  A  =  cost(S’)  -  cost(S). 

3.1.3  If  A  <  0  (downhill  move), 

Set  S  =  S’. 

3.1.4  If  A  >  0  (uphill  move), 

Set  S  =  S’  with  probability  e'A/T. 

3.2  Set  T  =  rT  (reduce  temperature). 

4.  Return  S. 

Source:  JOHNSON  ET  At.  [Ref  131. 


Figure  3.  Generic  Simulated  Annealing  Algorithm. 


•  there  is  limited  application  for  the  formulated  problem  (if  the  problem 
formulation  has  wide  applicability  efforts  toward  developing  an  effective 
problem  specific  algorithm  should  be  more  beneficial); 

•  the  solution  space  is  well  understood  (since  the  annealing  algorithm  is  highly 
parameter  dependent,  a  strong  understanding  of  the  solution  space  is  critical 
for  ensuring  convergence  of  the  algorithm  to  a  global  optimum). 

In  a  proposal  to  the  U.S.  Army  TRADOC  Analysis  Command  [Ref.  3],  Ingber 
proposed  that  a  variation  of  the  conventional  simulated  annealing  algorithm,  referred 
to  as  Very  Fast  Simulated  Reannealing,  could  be  used  effectively  to  fit  systems  of 
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differential  equations  to  battle  data  from  various  high  resolution  simulations  and 
NTC.  Our  initial  experiments  with  simulated  annealing  were  designed  to  use  a 
somewhat  simplified  version  of  Ingber’s  algorithm  to  assess  his  claim.  We  first 
attempted  to  fit  individual  battle  trajectories  and  various  subsets  of  battle  data  sets 
to  a  system  of  differential  equation  of  the  following  form  (our  difference  equation 
model  was  developed  later): 
dRTANK 


dt 

dRAPC 

dt 

dBTANK 


=  at 


rBTANK  +  a. 


,BAPC  +  a. 


vBTOW 


~  aRAj,cjtTANK  RTANK  +  ampcMPC  RAPC  +  oRAPC^TOW  BTOW 


£  aBTANK,RfANK  RTANK  +  a  BTANKJiAPC^^ 


dBAPC 

dt 

dBTOW 

dt 


-  aRAPCJtTANK  RTANK  +  oBAPCJiAPC  RAPC 


~  a BTOW, RTANK  RTANK  +  <^BrQWJtAPC  RAPC 


where  RTANK  denotes  the  number  of  red  tanks  involved  in  the  battle  at  time  t, 
a rtank. btank 's  ^  attrition  rate  of  red  tanks  by  blue  tanks,  and  similarly  for  the 
other  variables  and  coefficients  (note  here  that  APC  denotes  an  armored  personnel 
carrier,  and  TOW  denotes  the  Tube-Launched,  Optical-tracked,  Wire-guided  antitank 
weapon  system  of  the  U.S.  Army). 

After  numerous  experimental  runs  of  the  simulated  annealing  algorithm  it 
became  very  clear  that  regardless  of  the  parameters  used  in  the  algorithm,  when  a 
small  number  of  battles  were  used,  the  results  were  quite  unstable.  In  particular, 
starting  with  one  seed  for  the  pseudo-random  number  generator  we  would  get  one 
solution  and  starting  with  a  different  seed  we  would  arrive  at  a  different  solution 
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using  the  same  data.  When  a  larger  number  of  battles  were  used  the  computer  time 
required  became  unreasonable  (on  the  order  of  24  hours  on  a  386-based  PC  running 
at  33  Mhz  when  fitting  ten  battle  trajectories).  Based  on  these  results  we 
temporarily  abandoned  simulated  annealing  and  searched  for  quicker  and  more 
stable  techniques.  We  discuss  several  such  techniques  in  the  remainder  of  this 
chapter.  We  found  during  our  experiments  with  various  fitting  methods  that  it  was 
reasonable  to  specify  relative  fire  power  indices  for  each  weapon  system,  thus 
allowing  us  to  fit  a  simpler  system  of  differential  equations  given  below: 

2. 

dt 

dY  _  ,  v  7 

dt  y 

where, 

X  ~  f rtankRTANK  +  f RApcRAPC , 

Y  -  f btank^TANK  +  f^BAPC  +  fBmwBTOW. 

Experimental  runs  using  simulated  annealing  for  fitting  this  system  of  equations 
were  both  stable  and  much  quicker  than  those  experienced  previously.  However,  we 
observed  that  the  simulated  annealing  process  simply  converged  to  the  regression  fits 
discussed  in  Section  E,  and  took  about  10000  times  longer  to  do  so.  Thus,  simulated 
annealing  was  found  to  be  an  unsuitable  technique  for  fitting  systems  of  equations 
of  the  form  in  which  we  were  interested. 

C.  STEEPEST  DESCENT 

While  trying  to  determine  why  simulated  annealing  was  not  working  well,  we 
made  two  observations  regarding  our  model.  First,  Lanchester  models  assume  that 
every  weapon  system  on  the  battlefield  can  see  and  engage  every  other  system  on  the 
battlefield.  This  was  not  true  of  either  the  battles  observed  at  NTC  or  those 
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generated  by  Janus(A).  Second,  since  the  models  under  consideration  were  attrition 
models,  their  solution  spaces  were  convex,  and  thus  a  least  squares  measure  (i.e.  the 
sum  of  the  square  differences)  of  the  fit  of  the  data  to  these  models  would  also 
induce  a  convex  space.  Thus  a  simpler  fitting  technique  such  as  a  standard  steepest 
descent  algorithm  should  be  effective. 

At  this  point  we  consulted  with  analysts  at  the  TRADOC  Analysis  Command 
(TRAC)  and  determined  that  although  battle  data  showing  which  systems  could  see 
and  engage  which  other  systems  was  not  currently  available,  the  Janus(A)  system 
could  be  modified  to  produce  such  data.  Thus,  we  proceeded  by  generating 
simulated  battle  trajectories  from  mathematical  models  to  use  in  testing  our  fitting 
techniques  (see  Chapter  IV,  Section  A). 

We  began  testing  the  straight-forward  steepest  descent  algorithm  shown  in  Figure 
4,  using  this  simulated  data  (FORTRAN  code  is  in  Appendix  C). 

This  algorithm  proved  to  be  stable  with  respect  to  the  quality  of  solution,  but  not 
with  respect  to  time  of  solution.  That  is,  regardless  of  the  starting  solution  ,  the 
algorithm  would  converge  to  the  same  solution  for  a  given  accuracy  level;  however, 
the  time  required  to  do  so  varied  from  one  minute  up  to  four  hours.  Although  this 
algorithm  was  far  more  useful  than  the  simulated  annealing  approach,  the  idea  of 
using  non-linear  fitting  techniques  led  us  to  try  using  commercial  non-linear 
programming  software  to  fit  our  system  of  equations. 

D.  NON-LINEAR  OPTIMIZATION 

The  convex  nature  of  our  problem  encouraged  us  to  believe  that  we  might  be 
able  to  use  commercial  non-linear  solvers  to  fit  our  equations  to  the  battle  data.  We 
thus  formulated  our  problem  on  the  General  Algebraic  Modeling  System  (GAMS) 
using  the  MINOS  5.1  solver  (see  code  in  Appendix  D).  We  made  experimental  trials 
using  both  a  least  squares  measure  (i.e  Euclidean  distance)  and  a  Chebychev  (or 
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minimax)  criterion. 


Least  Squares  Objective  Function:  Min  V'  (6r2  +  6V2 

.  \  xi  t 

I 

Chebychev  Objective  Function:  Min  ^maxjl^j  +  |8y  |j 
Subject  to:  bx  k  (xi+1  -  xt)  -  (ayi  +  zx) 

by,  *  O’M  -  y)  -  &Xi  +  V 
Solve  for:  aopr  bopT>  (zx)opf  (zy)opT 


1.  Get  an  initial  solution  S. 

2.  Determine  an  initial  grid  definition  size. 

3.  While  not  yet  close  enough : 

3.1  Compute  cost  for  each  neighboring  point  on  the  grid. 

3.2  Set  S’  =  neighboring  solution  with  min  cost. 

3.3.  IF  IIS  -  S’ll  <  accuracy  THEN  refine  the  grid  size. 

4.  Return  S  =  S’. 


Figure  4.  Steepest  Descent  Algorithm. 


where  =  difference  between  the  predicted  attrition  and  the  actual  attrition. 

Our  results  were  encouraging.  The  solutions  were  stable  and  were  found 
consistently  in  just  under  one  minute  for  single  battles.  We  also  observed  graphically 
that  the  least  squares  criterion  provided  a  smaller  variance  on  the  fitted  coefficients. 

E.  LINEAR  REGRESSION 

Our  GAMS  formulation  with  a  least  squares  objective  function  was  essentially 
providing  a  simultaneous  linear  regression  fit  of  the  battle  trajectories  to  the 
equations  in  the  analytical  model  being  fit.  This  led  us  to  ask  when  this  would  be 
the  equivalent  or  nearly  equivalent  to  doing  linear  regression  on  each  of  the 
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individual  equations  separately.  The  answer  is  that  if  we  assume  the  coefficients  in 
the  equations  are  independent,  then  we  could  simply  fit  the  equations  separately 
using  regression. 

Considering  the  method  (Chapter  IV,  Section  A)  used  to  generate  our  simulated 
battle  trajectories  (which  held  the  attrition  and  noise  coefficients  constant  for  the 
duration  of  each  individual  battle),  this  independence  relationship  was  true  for  the 
individual  battle  trajectories  (this  would  not  be  the  case  for  the  data  in  general).  We 
did  regression  fits  on  the  individual  equations  and  found  the  two  methods  equivalent 
(i.e.  the  fitted  parameters  were  nearly  equal)  up  to  a  factor  that  could  be  explained 
by  noise  induced  by  the  generation  of  integer  data  points. 

Thus,  if  we  are  willing  to  make  the  assumption  that  the  equations  are 
independent  for  individual  battles,  we  reduce  our  fitting  time  from  about  one  minute 
per  battle  to  about  0.1  seconds  per  battle.  Whether  this  makes  sense  to  do  on  actual 
battle  trajectories  needs  to  be  examined  once  the  data  become  available. 

F.  CONCLUSION 

Our  experiments  with  the  non-linear  optimization  approach  of  using  a  GAMS 
model  with  a  least  squares  objective  function  indicate  that  this  approach  is  the  most 
effective  way  of  fitting  a  system  of  difference  equations  to  battle  trajectories, 
providing  consistent  solutions  in  a  timely  manner.  However,  experiments  should  be 
made  to  consider  whether  assuming  independence  between  the  coefficients  in  each 
equation  of  the  surrogate  model  is  justified,  thus  allowing  us  to  use  the  faster 
independent  regression  technique  discussed  in  the  last  section.  We  use  the 
independent  regression  technique  for  an  example  in  the  next  chapter. 
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IV.  EXAMPLE  OF  METHOD 


A.  THE  BATTLE  DATA  TRAJECTORIES 

To  test  the  methodology,  we  first  generated  battle  trajectory  data  based  on  a 
discrete  version  of  Lanchester’s  aimed  fire  model  with  additive  noise  term  (equation 
29  below).  Three  different  systems  of  equations  (all  of  this  same  form)  were  used 
for  this  purpose.  These  generating  systems  differed  only  in  the  fire  power  indices 
assigned  to  each  weapon  system.  Each  system  was  used  to  generate  two  separate 
collections  of  20  battles  each.  The  form  of  the  generating  equations  is: 

-  -ay.  + 

r..i  -  yu  -  ~BXn  +  (29) 

where 

Xn  =  f RTANK  FI  +  ffUPC  it  * 

Y.  -  /btank  I™*,  +  fsAPC  B^PCn  +fBTOWBTOWn  . 
and  A  and  B  are  constant  throughout  a  single  battle  trajectory.  Here,  fRTANK  is  the 
fire  power  index  for  the  red  tanks,  and  similarly  for  the  fire  power  indices  of  the 
remaining  weapon  systems  (values  shown  in  Table  2  below).  RTANK„  is  the  number 
of  red  tank  systems  active  in  the  battle  at  the  end  of  time  increment  n.  RTANKn  is 
computed  as  follows  from  Xn.  The  other  weapon  system  levels  (RAPC,  BTANK, 
BAPC,  and  BTOW)  are  all  computed  in  the  a  similar  manner.  At  the  end  of  each 
time  interval  [n,  n+  1],  the  attrition  sustained  by  the  X  force  is  divided  among  the 
two  X  weapon  systems  as  follows: 


RTANK  .  =  RTANK"  +  (X  . 

n  \  n+1 


tfpwrx  fRTANK 
(tfpwrx)(Xsys  -  1) 


RAPCntl  =  RAPCn  +  (Xn^  -  xn) 


-f 


RAPC 


(tfpwrx)  ( Xsys  -  1) 
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where 


tfPWTX  ~  f RTANK  +  f RAPC 

Xsys  s  the  number  of  different  weapon  system  in  X  force, 
and  similarly  for  the  Y  force.  We  then  round  RTANK^+j  and  RAPCn+1  down  to  the 
next  integer  value. 

Also,  [A  Zx  B  Zy]  is  distributed  multivariate  normal  with  variance-covariance 
matrix: 

0.00004  -0.00200  -0.00001  0.00010' 

-0.00200  0.50000  -0.00060  0.40000 
V  = 

-0.00001  -0.00060  0.00008  -0.00600 
0.00010  0.40000  -0.00600  1.00000 

and  mean  n  =  [  0.066  0.000  0.055  0.000  ]. 

Initial  force  levels  were 

RTANK(O)  =  80  BTANK(O)  =  60 

RAPC(0)  =  160  BAPC(0)  =  120 

BTOW(O)  =  40  , 


The  fire  power  indices  f  were  set  as  shown  in  Table  2: 


TABLE  2.  FIRE  POWER  INDICES. 


BATTLE  SET 

A  &  B 

C  &  D 

E  &  F 

RTANK 

0.90 

1.00 

0.90 

RAPC 

0.80 

0.90 

0.80 

BTANK 

1.00 

0.90 

1.00 

BAPC 

0.90 

0.80 

0.90 

BTOW 

0.80 

0.70 

0.25 
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A  realization  of  [A,  Zx,  B,  Zy]  was  chosen  and  a  battle  trajectory  consisting  of  12, 
five  minute  intervals  was  run  out  using  the  generating  equations  (29).  Another 
realization  was  then  chosen  and  another  battle  trajectory  was  generated.  In  this 
manner  20  battle  trajectories  were  generated  for  each  battle  data  set.  A  sample 
battle  trajectory  generated  with  these  equations  is  shown  in  Table  3. 


TABLE  3.  BATTLE  TRAJECTORY  SAMPLE. 


Battle  Trajectory  (Battle  1  Set  A) 

TIME 

RTANK 

RAPC 

BTANK 

BAPC 

BTOW 

0 

80 

160 

60 

120 

40 

5 

73 

152 

57 

117 

37 

10 

67 

145 

54 

114 

33 

15 

61 

138 

51 

111 

30 

20 

55 

132 

49 

108 

27 

25 

49 

126 

46 

106 

25 

30 

44 

120 

44 

103 

22 

35 

39 

114 

42 

101 

20 

40 

34 

109 

40 

99 

18 

45 

30 

103 

39 

97 

16 

50 

25 

98 

37 

96 

14 

55 

21 

94 

35 

94 

13 

60 

17 

89 

34 

93 

11 

B.  SELECT  MATHEMATICAL  SURROGATE  SYSTEM 

We  had  decided  to  use  the  Lanchester  aimed  fire  model  in  our  initial  tests. 
However,  we  had  to  decide  whether  to  include  the  additive  (random)  noise  term. 
The  concern  was  to  choose  the  model  that  would  best  represent  the  underlying 
structure  of  the  data  represented  by  the  generating  system.  The  best  candidate  was 
determined  by  running  a  small  sample  test.  Five  battles  were  selected  at  random 
from  battle  set  A  (see  Table  3)  and  were  fit  to  the  discrete  form  of  the  Lanchester 
aimed  fire  model  with  and  without  the  constant  noise  terms.  The  optimization 
model  (procedure)  used  for  estimating  the  parameters  in  the  following  surrogate 
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model  without  additive  noise  terms 

- *.  *  [-« * 

r..,  -  r,  -[-»  *  c,(r.)]X„,  (30) 

is  given  below. 

Proportional  Procedure  (no  additive  noise  term): 

Minimize:  £(6X2  +  6y2) 

i 

subject  to:  bx  ;>  (XM  -  X)  -  ( a  Y) 

\  *  VM  -  Y)  ~  0>X,) 

Solve  for:  dOFT>  ^  opt 
where  and  5y  are  as  in  Chapter  II. 

Additionally,  the  surrogate  model  with  the  additive  noise  term  (equation  27)  was 
fit  using  a  one-stage  and  two-stage  fitting  procedure.  The  one-stage  procedure  was 
the  non-linear  programming  technique  describe  in  Chapter  III. 

One-stage  Procedure: 

Minimize:  £  (6x2  + 

subject  to:  bx  z  (Xm1  -  X)  -  (a  Yf  +  zx) 

6>,  *  ~  Y)  ~  (bX,  + 

Solve  for:  dopv  bOPV  iXopJ>  z  >orr 

The  two-stage  procedure  used  the  same  basic  technique;  however,  in  stage  one  the 
noise  terms  were  fixed  at  zero  and  the  attrition  coefficients  were  fit  (this  stage  is  the 


39 


same  as  the  proportional  procedure).  In  the  second  stage  the  attrition  coefficients 
were  fixed  at  the  optimal  value  found  in  stage  one  and  the  noise  term  was  fit. 

Two-stage  procedure: 

Stage  One  Problem: 

Minimize: 
subject  to: 

Solve  for: 

Stage  Two  Problem: 

Minimize:  £(6*,2  +  V) 

I 

subject  to:  bx  z  (Xitl  -  X ) 

»„ 2  o-..,  - 1',) 

Solve /or.  z,oFf  Z)ofi 

The  results  are  shown  in  Tables  4,  5,  and  6,  where  cost  is  taken  to  be  the  square  root 
of  sum  of  square  differences  of  the  predicted  and  observed  attrition  from  the  model 
(equation  29). 


opt  Y,  +  Zz) 
{b  opr  Xi  +  Zy) 


Z(6,2  +  M 

\  >  (XM  -  X)  -  (ft  Y') 
\  *  Oj* i  -  Y)  -  (bX) 

&OPV  b  OPT 
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TABLE  4.  ONE-STAGE  FITTING. 


TABLE  5.  TWO-STAGE  FITTING. 


Results  of  Two-Stage  Fitting  Procedure 

BATTLE 

A-l 

A-2 

A-3 

> 

i 

-p. 

A-5 

COST 

4.675 

8.742 

3.720 

5.700 

7.938 

a 

0.061 

0.048 

0.067 

7 

0.030 

0.021 

0.056 

*-x 

0.042 

0.048 

0.037 

0.049 

b 

*, 

-3.519 

0.176 

-1.764 

-0.539 

3.082 

TABLE  6.  PROPORTIONAL  FITTING. 


Results  of  Proportional  Fitting  Procedure 

BATTLE 

A-l 

A-2 

A-3 

i 

< 

A-5 

COST 

5.378 

7.047 

4.352 

5.147 

10.991 

a 

1 

0.048 

■  : 

b 

; 

0.037 

Examining  the  costs  (sum  of  squares  error)  in  the  tables  above  we  concluded  that 
the  one-stage  fitting  procedures  was  the  best  (minimum  error)  of  the  three  examined. 
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We  were  however  interested  in  obtaining  a  fitted  model  that  best  represents  the 
underlying  structure  of  the  data  that  we  were  fitting.  Thus,  we  also  looked  at  the 
Euclidean-norm  of  the  difference  between  the  fitted  coefficients  and  those  used  to 
generate  the  data.  Let  a  =  [A,  Zx,  B,  Zy];  define  a  to  be  the  parameters  estimated 
by  the  fitting  procedure  (take  Zx  =  Zy  =  0  in  the  proportional  procedure).  Let 
\xz  ,  \lb,  \iz  ],  the  mean  values  of  the  parameters  used  in  the  generating 

equations  to  simulate  the  battle  data  sets.  We  then  compute 
I  i  -  I  -  /(<i  -  *  (ix  -  li2_f  *  {(  -  »,f  *  (i,  -  values  are 


shown  in  Table  7. 

TABLE  7.  NORM  OF  DIFFERENCES  IN  COEFFICIENT  VECTORS. 


Euclidean-norm  of  difference  in  Coefficients 

BATTLE 

A-l 

jjfjjj 

A-3 

AVG 

2-Stage 

2.59 

0.04 

1.96 

1.78 

3.10 

1.85 

1 -Stage 

2.24 

1.88 

1.56 

2.22 

2.41 

2.06 

Pro. 

0.93 

0.25 

0.21 

1.45 

0.58 

0.63 

We  also  examined  the  battle  trajectories  for  each  of  these  five  fitted  surrogate 
models  graphically,  by  plotting  the  force  levels  versus  time  for  the  fitted  system  and 
the  generating  system  (see  Figures  5-7  for  an  example). 

As  a  result  of  the  analyses  we  concluded  that  although  the  one-step  fit  was 
providing  the  best  fit  in  the  least  squares  sense;  when  we  examined  the  results  using 
the  proportional  fit  on  each  individual  equation  (in  equation  28)  independently 
(Figure  8)  the  fit  to  the  underlying  structure  (i.e.  the  generating  equations)  was 
providing  a  graphical  representation  almost  as  good  as  the  one-stage  procedure. 
Since  the  proportional  procedure  used  was  much  quicker  than  the  one-stage 
procedure  we  decided  to  use  it.  The  following  optimization  model  is  the  procedure 
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Time  Sequence  Pi 


TRGEN  RED 
TR1SG  RED 


Figures.  Battle  Trajectory  Plots.  (1-Stage  vs  Gen) 


Figure  6.  Battle  Trajectory  Plots.  (2-Stage  vs.  Gen) 
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used. 

Proportional  independent  procedure  (PRO-I): 


a  = 


i-0 


y. 


i  / 


N 


n-i  /  y  _  y  \ 

[  ■'i-l 


b  = 


j  =0 


x.. 


«  ) 


N 


We  therefore  used  the  discrete  difference  equation  model  without  additive  noise 
(equation  30)  for  the  remainder  of  the  analysis. 

C.  FITTING  THE  MODEL  TO  THE  DATA 

Based  on  the  analysis  in  the  preceding  section,  the  model  was  fit  to  the  battle 
data  sets  obtained  from  the  generating  systems  using  the  proportional  independent 
method.  The  coefficients  fit  to  the  first  five  replications  of  each  data  set  are  shown 
in  Table  8  below. 
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TABLE  8.  FITTED  COEFFICIENTS. 


Sample  of  Fitted  Coefficients 

Battle 

Set  A 

Set  B 

Set  C 

Set  D 

Set  E 

Set  F 

H 

IB 

Ull 

0.0493 

0.0499 

0.0406 

0.0582 

0.0499 

0.0417 

0.0414 

0.0466 

2:  & 

b 

0.0489 

0.0530 

0.0383 

0.0548 

■ 

msm 

B 1 

3:  d 

b 

0.0533 

0.0396 

0.0409 

0.0515 

1 

■ 

M 

4:  a 

b 

1 

0.0573 

0.0499 

R 

0.0479 

0.0464 

5:  d 

b 

■ 

0.0572 

0.0470 

— 

R 

■ 

0.0476 

0.0446 

Mean: 

a 

b 

RjHK  1 

Ki§ 

1 

I 

0.0471 

0.0555 

0.0462 

0.0448 

0.0479 

0.0452 

D.  EMPIRICAL  DISTRIBUTION  FUNCTIONS 

Graphical  analysis  of  the  fitted  coefficients  from  the  battle  data  sets  (of  20  battle 
trajectories  each)  show  no  clear  fit  to  a  known  probability  distribution,  thus  we 
turned  to  the  empirical  distributions.  The  empirical  cumulative  distribution  functions 

for  the  ratio  of  the  estimated  coefficients  —  where  determined  for  each  battle  (see 

a 

Figure  9  below). 

Using  the  procedure  outlined  in  Chapter  II  we  generated  the  probability  surfaces 
numerically  over  square  intervals  and  calculated  the  difference  between  surface 
heights  at  the  corner  of  each  square  interval  (X  <  Xo  <  X+l,  Y  <y0  <  Y+l}.  The 
observed  values  of  our  test  statistic  T  (volume  between  surfaces  compared)  of  the 


48 


battle  set  comparisons  are  listed  in  Table  9.  The  notched  box  plots  in  Figure  10 
provide  a  graphical  portrayal  of  samples  of  size  50  of  the  randomization  distribution 
of  these  differences  and  the  location  of  the  observed  test  statistic  within  this 
distribution.  Examination  of  these  box  plots  clearly  delineates  the  differences 
between  the  dissimilar  battle  data  sets,  and  the  similarity  of  the  equivalent  battle 
data  sets  (i.e.  A  and  B,  C  and  D,  and  E  and  F  as  defined  in  Table  3). 


TABLE  9.  OBSERVED  VALUE  OF  THE  TEST  STATISTIC  T. 


Observed  Value  of  Test  Statistic 

vs 

B 

C 

D 

E 

F 

A 

201 

2940 

2849 

1079 

928 

B 

3023 

2932 

1162 

1011 

C 

443 

1860 

2012 

D 

1769 

1927 

E 

283 

E.  ANALYTICAL  COMPARISON  OF  BATTLE  DATA  SETS 


Having  calculated  the  differences  between  probability  surfaces  generated  from 
the  empirical  distribution  functions  of  random  samples  of  the  ratios  of  fitted 
coefficients,  we  now  use  the  50  samples  obtained  from  the  randomization  distribution 
of  the  surfaces  differences  and  the  observed  value  of  our  test  statistic  (T)  to  compare 
the  battle  data  sets  from  the  two  simulation  models  being  compared.  To  do  so  we 
simply  observe  the  relative  location  of  the  observed  statistic  in  the  randomization 
distribution  by  calculating  the  percentage  of  realization  values  greater  then  the 
observed  value.  This  is  equivalent  to  determining  the  percentile  in  which  it  lies.  If 
the  percentage  is  large  then  we  do  not  reject  a  hypothesis  that  the  two  battle  data 
sets  came  from  similar  (possibly  the  same)  simulation  models.  However,  if  the 
percentage  is  small  then  we  reject  that  hypothesis.  Since  our  trials  resulted  in 
percentages  well  out  in  the  tails  of  the  randomization  distributions  (see  Table  10), 


50 


we  did  not  postulate  any  conclusions  with  respect  to  an  exact  cut  off  percentage  for 
rejection. 


TABLE  10.  ANALYTICAL  COMPARISONS  OF  RESULTS  (THE  NUMBERS 
INDICATE  THE  PERCENTAGE  OF  SAMPLE  VALUES  OF  THE 
RANDOMIZATION  DISTRIBUTION  GREATER  THAN  THE  TEST 
STATISTIC). 


Analytical  Comparison 

B 

C 

D 

E 

F 

DO  NOT 

REJECT 

REJECT 

REJECT 

REJECT 

REJECT 

0.96 

0.00 

0.00 

0.26 

0.14 

REJECT 

REJECT 

REJECT 

REJECT 

KB  J 

0.00 

0.00 

0.08 

0.08 

DO  NOT 

REJECT 

REJECT 

KC,J 

REJECT 

0.74 

0.10 

0.04 

REJECT 

REJECT 

KD,J 

0.04 

0.0.'' 

DO  NOT 

KE,_) 

REJECT 

0.90 
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V.  CONCLUSIONS  AND  RECOMMENDATIONS 

A.  INTRODUCTION 

This  thesis  has  presented  a  method  for  comparing  highly  complex  combat  models 
using  a  simple  analytical  surrogate  model.  In  doing  so,  we  have  evaluated  a  number 
of  fitting  techniques  and  other  comparison  techniques,  and  we  have  developed  a 
discrete  analog  to  Lanchester’s  aimed  fire  model. 

B.  SURROGATE  MODELS 

We  have  concentrated  here  on  Lanchestei-like  models  of  combat;  however,  other 
models  could  be  studied.  Non-homogeneous  models  which  include  area  fire  weapons 
(such  as  artillery  and  mortars)  were  not  studied  in  this  thesis  but  should  be  studied 
if  this  methodology  is  to  be  useful  as  an  analytical  tool.  Models  specifically 
addressing  human  factors  should  also  be  explored. 

C.  FITTING  TECHNIQUES 

Although  we  found  that  the  simultaneous  least  squares  (GAMS  model)  technique 
most  effectively  represented  the  underlying  structure  of  battle  data  sets  fitted  to  our 
model,  all  of  the  fitting  techniques  considered  have  utility  for  certain  types  of 
problems.  Simulated  annealing  is  a  very  time  consuming  procedure,  but  could  be 
useful  if  the  simulation  models  being  studied  allow  reinforcement.  A  useful 
algorithm  for  fitting  non-linear,  non-convex  continuous  functions  using  simulated 
annealing  is  discussed  in  Brooks  and  Verdini  [Ref.  15]. 

GAMS  provided  a  very  flexible  experimental  tool  for  fitting  our  models. 
Numerous  possible  fitting  techniques  can  be  written  in  GAMS  in  an  efficient  manner. 

We  expect  that  attempts  to  use  multiple  regression  (as  in  our  proportional 
independent  procedure)  with  more  complex  models  will  experience  difficulty  with  the 
independence  assumption.  It  has,  however,  provided  a  powerful  tool  for  fitting  a 
large  number  of  battles  in  a  very  short  time. 
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D.  RECOMMENDATION 

Comparing  high  resolution  models  of  combat  continues  to  be  a  challenging 
problem.  We  have  seen  here  that  it  is  feasible  to  use  very  simple  analytic  surrogate 
models  to  compare  machine  generated  battle  data  sets.  There  is  also  a  lot  of  work 
to  be  done  in  applying  this  methodology  to  actual  data  sets  generated  by  Janus(A) 
and  at  NTC,  and  in  investigating  various  surrogate  models.  Successful  comparison 
of  these  two  highly  complex  models  of  combat  should  be  of  great  benefit  to  the  Army 
analysis  community. 
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APPENDIX  A 

MATRIX  SOLUTION  TO  DISCRETE  MODEL 

In  Chapter  III  we  developed  the  2x2  case  of  the  discrete  difference  equation 
model.  In  doing  so  we  combined  the  fire  power  from  the  various  weapon  systems 
into  one  term  by  using  the  fire  power  indices.  Suppose  instead  we  wish  to  consider 
a  model  where  the  attrition  of  each  weapon  system  and  its  contribution  to  the 
attrition  of  all  other  weapon  systems  are  considered  separately.  The  discrete  model 
for  this  system  would  be  given  by: 

Xj.k+i  -  X,k  =  -an  Ylk  -a12  Y2k  •  •  •  -alm  Ym  k , 

^2,k  + 1  *  X2.k  =  -a21  Ylk  -a22  Y2  k  •  •  •  -a2m  Ym  k  , 

• 

Xn,k+i  '  Xn  k  —  -anl  Y[  k  -an2  Y2k  •  •  •  -anm  Ymk , 

^i.k+i  ‘  Ylik  =  -bn  Xl  k  -bI2  X2k  •  •  •  -bln  X,^  , 

^2,k  +  i  *  Y2k  =  -b21  XI  k  -b22  X2k  •  •  •  -b2n  Xnk  , 

^m.k  +  l  *  ^m.k  “  ”bm]  Xj  k  "bm2  X2  k  •  •  •  *bmn  Xn  k  . 
where, 

X,  k  =  the  number  of  X  force  weapon  type  i,  active  in  the 
battle  at  the  end  of  time  increment  k, 

Yj  k  s  the  number  of  Y  force  weapon  type  j,  active  in  the 
battle  at  the  end  of  time  increment  k, 

,  3  the  rate  at  which  Y  force  weapon  type  j  kills 
X  force  weapon  type  i, 

b, j  s  the  rate  at  which  X  force  weapon  type  i  kills 
Y  force  weapon  type  j. 
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or  in  matrix  notation: 


*u*l' 

1  0  0  •  •  •  0  -an  •  •  •  -fllM 

V 

*2*1 

0  1  0  •  •  •  0  -al2  •  •  •  -a^ 

* 

V. 

0  •  •  ‘  01  - anl  •  •  • 

X, 

-bn  -  •  •  -bln  1  0  0  •  •  •  0 

-b2\  •  •  •  -b2n  0  1  0  •  •  •  0 

.  .  ^ 
V 

A 1  •••  -Kn  0  •••  0  1 

^  • 

> 

In  block  notation  we  can  write: 


I 

A 

X 

V 

B 

I 

Yk 

JM 


R  Fk  ;  where 


R  = 


/  A 
B  I 


The  solution  to  this  system  is: 


F(k)  =  Rk  F(0) . 

Associated  with  R  are  its  eigenvalues  {A,,  A2,  •  •  •  ,  An  +  m}  and  their  corresponding 
eigenvectors  {A,,  A2,  •  •  •  ,  An  +  m}.  These  eigenvectors  form  a  basis  for  Rn  +  m  so  we 
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can  then  write  F(0)  as  a  linear  combination  of  these  eigenvectors: 

F(  0)  =  ci  Aj  +  c2  Aj  +  •  •  •  +  cn»m  A„+m. 

So  finally,  the  solution  to  the  discrete  system  can  be  written: 

TO  =  ci(Xi)*Ai  +  *  *  *  '  +c»*m  (*»*«)*  A» 
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APPENDIX  B 

SIMULATED  ANNEALING  CODE 

PROGRAM  siman 

C**  Simulated  Annealing  Algorithm  implemented  for  the 
C**  discrete  analog  to  Lanchester’s  Aimed  Fire  Model. 

REAL  cred(  100),cblue(  100),red(  100),blue(  100), c(4), 
&  acost,ccost,fnacc,fniter, pace, temp, cool, minpct 

INTEGER  n, max  eye 
COMMON/vars/cred,cblue, red, blue, n 

INTEGER  nacc, niter 

REAL  tcost,mesh(4),a(4),at(4) 

DOUBLE  PRECISION  dseed 

dseed  =  1 8564245. OdO 

DO  31  i  =  1,4 
a(i)  =  0.0 
c(i)  =  0.0 
31  CONTINUE 

mesh(l)  =  0.0001 
mesh(2)  =  0.001 
mesh(3)  =  0.0001 
mesh(4)  =  0.001 

temp  =  100.0 
cool  =  0.95 
max  eye  =  100 
minpct  =  0.01 

CALL  readit 
CALL  evalfn(a,acost) 
ccost  =  acost 
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C“  Begin  Simulated  Annealing  Algorithm 
DO  36  k  =  1  ,  max  eye 
nacc  =  0 
niter  =  5000 
DO  35  j  =  1  ,  niter 
CALL  nbhd(a,at,mesh,dseed) 

CALL  evalfn(at,tcost) 

IF  (tcost  .le.  acost)  THEN 
DO  33  i  =  1  ,  4 
a(i)  =  at(i) 

IF  (tcost  .It.  ccost)  c(i)  =  a(i) 

33  CONTINUE 
acost  =  tcost 

IF  (acost  .It.  ccost)  ccost  =  acost 
nacc  =  nacc  +  1 
ELSE 

CALL  unif(dseed,u) 

IF  (u  .gt.  exp(-temp*(tcost-acost)))  THEN 
DO  34  i  =  1  ,  4 
a(i)  =  at(i) 

34  CONTINUE 
acost  =  tcost 
nacc  =  nacc  +  1 

ENDIF 

ENDIr 

35  CONTINUE 

fnacc  =  FLO  AT(  nacc) 
fniter  -  FLOAT(niter) 
pace  =  fnacc  /  fniter 
C“  Geometric  Cooling 
temp  =  cool  *  temp 
IF  (pace  .It.  minpet)  THEN 

PRINT  *,’End  conditions  met  in  \k,’  cycles’ 
STOP 
ENDIF 

36  CONTINUE 

PRINT  ‘/Maximum  iterations  reached’ 

201  FORMAT(lx,’A:\4(F12.5),/,lx,’C:\4(Fl2.5)) 

202  FORM.  \T(lx,’ACOST:\F  12.5,’  CCOST:\F12.5, 
&'  PCT  ACCEPTED:',F8.3,/) 

END 
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SUBROUTINE  readit 


REAL  cred(100),cblue(100),red(100),blue(100) 
INTEGER  n 

COMMON/vars/cred,cblue, red, blue, n 

INTEGER  ulist,ufpwr,udata,uout,d(5),numobs,i,time 
REAL  fpwr(5) 

CHARACTER *9  bfile 

ulist  =  8 

ufpwr  =  9 

udata  =  10 

uout  =  11 

OPEN(UNIT = ulist,FILE  =  ’slist.dat’, STATUS  =  ’old’) 
OPEN(UNIT  =  ufpwr, FILE  =  ’fpwr.dat’, STATUS  =  ’old’) 
DO  1  i  =  1,5 
READ(  ufpwr,  101)  fpwr(i) 

1  CONTINUE 
n  =  1 

2  CONTINUE 

C**  Read  list  of  Battle  Data  Set  Files  to  be  used. 
READ(ulist,  102,END  =  4)  bfile 
OPEN(UNIT  =  udata, FILE  =  bfile,STATUS  =  ’old’) 
READ(udata,103)  numobs 
DO  3  i  =  n  ,  n  +  numobs-1 

READ(udata,104)  time,d( l),d(2),d(3),d(4),d(5) 
red(i)  =  FLOAT(d(l))*fpwr(l)  + 
FLOAT(d(2))*fpwr(2) 

blue(i)  =  FLO  AT(d(3))  *  fpwr(3 )  + 
FLOAT(d(4))*fpwr(4)  + 

&  FLOAT(d(5))*fpwr(5) 

IF  (i  .gt  .  1)  THEN 
cred(i-l)  =  red(i)  -  red(i-l) 
cblue(i-l)  =  blue(i)  -  blue(i-l) 

ENDIF 

3  CONTINUE 

n  =  n  +  numobs-1 
GOTO  2 

4  CONTINUE 
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OPEN(  UNIT  =  uout,FILE  =  ’sa.out’,  STATUS  =  ’unknown’) 
DO  5  i  =  1  ,  n-1 

WRITE(uout,  105)  i,red(i),cred(i),blue(i),cblue(i) 

5  CONTINUE 
n  =  n-1 

101  FORMAT(F10.7) 

102  FORMAT(A9) 

103  FORMAT(I5) 

104  FORMAT(6I5) 

105  FORMAT(I3,5(2x,F10.7)) 

RETURN 

END 

SUBROUTINE  evalfn(c,ecost) 

C**  Sum  of  Square  Differences  (Linear  Lanchester) 

REAL  cred(100),cblue(100),red(100),blue(100) 
INTEGER  n 

COMMON/vars/cred,cblue,red,blue,n 

REAL  ecost,r,b,c(4) 

ecost  =  0.0 
DO  201  i  =  1  ,  n 

r  =  (cred(i)  -  (c(l)*blue(i)  +  c(2)))**2 
b  =  (cblue(i)  -  (c(3)*red(i)  +  c(4)))**2 
ecost  =  ecost  +  r  +  b 
201  CONTINUE 

ecost  =  SQRT(ecost) 

RETURN 

END 
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SUBROUTINE  nbhd(b,bt,m, dseed) 

C**  Chooses  random  neighboring  solution. 

REAL  b(4),bt(4),m(4),u,s 

DOUBLE  PRECISION  dseed 


CALL  unif(dseed,u) 

CALL  unif(dseed,s) 
j  =  INT((u*4.0)  +  0.5) 
sign  =  (s-0.5)/abs(s-0.5) 
bt(j)  =  b(j)  +  sign*m(j) 

RETURN 

END 

SUBROUTINE  UNIF(DSEED,U) 

C**  Uniform  (0,1)  Pseudo  Random  Number  Generator  (Lewis  & 
C**  Orav  [Ref.  16]) 

INTEGER  I 
REAL  U 

DOUBLE  PRECISION  DENOM, DSEED 

DATA  DENOM/2147483647.DO/ 

DSEED  =  DMOD(  16807.D0*DSEED, DENOM) 

U  =  DSEED/DENOM 

RETURN 

END 
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APPENDIX  C 

STEEPEST  DESCENT  CODE 


PROGRAM  steepd 

C**  Steepest  Descent  Algorithm 

REAL  cred(  100),cblue(  100),red(  100),blue(  1 00),c(4), 

&  acost,ccost,fnacc,fniter, pace, temp, cool, minpct 

INTEGER  n,  max  eye 
COMMON/vars/cred,cblue,red,blue,n 

INTEGER  nacc, niter, tmax, term 

REAL  tcost,mesh(4),a(4),at(4),amin(4),amax(4) 

DO  31  i  =  1  ,  4 
a(i)  =  0.0 
c(i)  =  0.0 
31  CONTINUE 

mesh(l)  =  0.0001 
mesh(2)  =  0.001 
mesh(3)  =  0.0001 
mesh(4)  =  0.001 
amin(l)  =  -1.0 
amax(l)  =  1.0 
amin(2)  =  -5.0 
amax(2)  =  5.0 
amin(3)  =  -1.0 
amax(3)  =  1.0 
amin(4)  =  -5.0 
amax(4)  =  5.0 

precis  =  0.0001 
refine  =  0.60 
max  eye  =  20 
maxit  =  2500 


CALL  readit 
CALL  evalfn(a,acost) 
pcost  =  acost 
ncyc  =  1 
it  =  1 


C  Top  of  while  loop . 

C**  Finds  Direction  of  greatest  improvement  in  the  Cost  Fen. 

1000  delmax  =  O.OdO 
DO  1001  term  =  1,4 

a(term)  =  a(term)  +  mesh(term) 

IF  (a(term)  .le.  amax(term))  THEN 
CALL  evalfn(a,tcost) 
delta  =  acost  -  tcost 
If  (delta  .gt.  delmax)  THEN 
delmax  =  delta 
tmax  =  term 
sign  =  l.OdO 
ENDIF 
ENDIF 

a(term)  =  a(term)  -  2.0d0*mesh(term) 

IF  (a(term)  .ge.  amin(term))  THEN 
CALL  evalfn(a, tcost) 
delta  =  acost  -  tcost 
IF  (delta  .gt.  delmax)  THEN 
delmax  =  delta 
tmax  =  term 
sign  =  -l.OdO 
ENDIF 
ENDIF 

a(term)  =  a(term)  +  mesh(term) 

1001  CONTINUE 

C**  If  an  improving  direction  is  found  move  in  that  Direction 
C**  on  grid  [mesh]  square. 

IF  ((delmax  .gt.  0.0)  .and.  (it  .le.  maxit))  THEN 
a(tmax)  =  a(tmax)  +  sign*mesh(tmax) 

CALL  evalfn(a, acost) 
it  =  it  +  1 
GOTO  1000 

C**  If  change  in  cost  during  this  cycle  less  then  set  Precision 
C**  We  are  done! 
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ELSEIF  ((pcost-acost)  .It.  precis)  THEN 
GOTO  1003 

C**  If  num  of  cycles  is  <  maximum  cycle  limit  =  >  Refine  Grid 
.It.  max  eye)  THEN 

ncyc  =  ncyc  +1 
pcost  =  acost 
DO  1002  term  =  1,  4 
mesh(term)  =  mesh(term)  *  refine 

1002  CONTINUE 
it  =  1 

GOTO  1000 
ENDIF 

1003  CONTINUE 

Print  ’ncyc, acost,  ncyc, acost, it, maxit 
Print  *,’a’,a 
PRINT  Vend’ 

WRITE(*,392)  ’  COST, FLO AT(acost) 

390  FORMAT(lx,A6,8x,A12) 

391  FORMAT(lx,A25,2x,Fl0.7) 

392  FORMAT(lx,A5,lQx,F12.7) 

36  CONTINUE 
END 


SUBROUTINE  readit 

REAL  cred(100),cblue(100),red(100),blue(100) 
INTEGER  n 

COMMON/vars/cred,cblue, red, blue, n 

INTEGER  ulist,ufpwr,udata,uout,d(5),numobs,i,time 
REAL  fpwr(5) 

CHARACTER*9  bfile 

ulist  =  8 
ufpwr  =  9 
udata  =  10 
uout  =  1 1 

OPEN(UNIT  =  ulist, FILE  =  ’slist.dat’,  STATUS  =  ’old’) 
OPEN(UNIT  =  ufpwr.FILE  =  ’fpwr.dat’, STATUS  =  ’old’) 


ELSEIF  (ncyc 
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DO  1  i  =  1  ,  5 
READ(ufpwr,101)  fpwr(i) 

1  CONTINUE 
n  =  1 

2  CONTINUE 
READ(ulist,102,END =4)  bfile 
OPEN(UNIT = udata,FILE  =  bfile, STATUS  =  ’old’) 

READ(udata,103)  numobs 
DO  3  i  =  n  ,  n  +  numobs-1 
READ(udata,104)  time,d(l),d(2),d(3),d(4),d(5) 
red(i)  =  FLOAT(d(l))*fpwr(l)  + 
FLOAT(d(2))*fpwr(2) 

blue(i)  =  FLOAT(d(3))*fpwr(3)  + 
FLOAT(d(4))*fpwr(4)  + 

&  FLOAT(d(5))*fpwr(5) 

IF  (i  .gt  .  1)  THEN 
cred(i-l)  =  red(i)  -  red(i-l) 
cblue(i-l)  =  blue(i)  -  blue(i-l) 

ENDIF 

3  CONTINUE 

n  =  n  +  numobs-1 
GOTO  2 

4  CONTINUE 

OPEN(UNIT  =  uout,FILE  =  ’sa.out’, STATUS  =  ’unknown’) 
DO  5  i  =  1  ,  n-1 

WRITE(uout,105)  i,red(i),cred(i),blue(i),cblue(i) 

5  CONTINUE 
n  =  n-1 

101  FORMAT(F10.7) 

102  FORMAT(A9) 

103  FORMAT(I5) 

104  FORMAT(6I5) 

105  FORMAT(I3,5(2x,F10.7)) 

RETURN 

END 
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SUBROUTINE  evalfn(c,ecost) 

C**  Sum  of  Square  Differences  (Linear  Lanchester  Model) 

REAL  cred(  100),cblue(  100),red(  100),blue(  100) 
INTEGER  n 

COMMON/vars/cred,cblue, red, blue, n 

REAL  ecost,r,b,c(4) 

ecost  =  0.0 
DO  201  i  =  1  ,  n 

r  =  (cred(i)  -  (c(l)’blue(i)  +  c(2)))**2 
b  =  (cblue(i)  -  (c(3)*red(i)  +  c(4)))**2 
ecost  =  ecost  +  r  +  b 
201  CONTINUE 

ecost  =  SQRT(ecost) 

RETURN 

END 
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APPENDIX  D 

GAMS  REGRESSION  MODEL 


STITLE  Least  Squares  Fit  of  Battle  Data  to  Analytical  Model 
SSTITLE  (Linear  Lanchester  Form  -  2-stage  fit) 

* . GAMS  OPTIONS  AND  DOLLAR  CONTROL  OPTIONS 


SOFFUPPER  OFFSYMXREF  OFFSYMLIST 


OPTIONS  LIMCOL  =  0,  LIMROW  =  0,  SOLPRINT  =  OFF 
OPTIONS  RESLIM  =  50,  ITERLIM  =  10000,  OPTCR  =  0.0; 
. DEFINITIONS  AND  DATA - - 


SETS 

I  increments  /l*  12/ 


W  weapons  / 

RTANK 
RAPC 
BTANK 
BAPC 
BTOW  /; 


PARAMETER  FP(W)  fire  power  / 

RTANK  .920 
RAPC  .780 
BTANK  1.000 
BAPC  .850 
BTOW  .750  /; 
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TABLE  CHANGE(I.W)  changes  in  wpn  level  in  inc 


RTANK  RAPC  BTANK  BAPC  BTOW 

1  7  8  4  4  4 

2  8  8  3  4  4 

3  6  8  4  3  4 

4  6  7  3  3  3 

5  6  6  2  3  3 

6  6  6  3  3  3 

7  5  6  2  2  2 

8  5  6  2  2  3 

9  4  5  2  2  2 

10  5  5  2  2  2 

11  4  5  12  1 

12  4  4  1  1  2; 

TABLE  LEVEL(I,W)  level  of  wpn  at  start  of  inc 


RTANK  RAPC  BTANK  BAPC  BTOW 


1 

80 

160 

60 

120 

40 

2 

73 

152 

56 

116 

36 

3 

65 

144 

53 

112 

32 

4 

59 

136 

49 

109 

28 

5 

53 

129 

46 

106 

25 

6 

47 

123 

44 

103 

22 

7 

41 

117 

41 

100 

19 

8 

36 

111 

39 

98 

17 

9 

31 

105 

37 

96 

14 

’0 

27 

100 

35 

94 

12 

11 

22 

95 

33 

92 

10 

12 

18 

90 

32 

90 

9  ; 

PARAMETER  RED(I)  red  firepower  at  increment  i; 


RED(I)  =  FP( "RTANK")  *  LEVEL(I,"RTANK") 

+  FP("RAPC")  *  LEVEL(I,"RAPC"); 

PARAMETER  CRED(I)  change  in  red  firepower  in  increment; 

CRED(I)  =  FPC'RTANK")  *  CHANGE(I,"RTANK") 

+  FPC'RAPC")  •  CHANGE(I,"RAPC"); 
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PARAMETER  BLUE(I)  blue  firepower  at  increment  i; 

BLUE(I)  =  FP('BTANK')  *  LEVEL(I,"BTANK") 
+  FP('BAPC')  *  LEVEUI.'BAPC") 

+  FPfBTOW")  *  LEVEL(I,"BTOW”); 


PARAMETER  CBLUE(I)  change  in  firepower  of  blue  in  increment; 

CBLUE(I)=  FP("BTANK")  *  CHANGE(I,"BTANK") 

+  FP("BAPC")  *  CHANGE(I,"BAPC") 

+  FP("BTOW")  *  CHANGEa,"BTOW"); 


. MODEL . 

POSITIVE  VARIABLES 

A  coef  on  blue  level  in  equation  red 

C  coef  on  red  level  in  equation  blue 

BL(I)  residual  in  blue  equation 

RD(I)  residual  in  red  equation; 


VARIABLES 

B  additive  noise  term 

D  additive  noise  term; 


VARIABLE 

COST  sum  of  squared  error; 


EQUATIONS 

OBJ 

DREDl(I) 

DBLUEl(I) 

DRED2(I) 

DBLUE2(I) 


minimize  sum  of  squared  error 
sq  error  increment  no  noise  term 
sq  error  increment  no  noise  term 
sq  error  in  increment  with  noise  term 
sq  error  in  increment  with  noise  term; 


>>>>>>>>>>  MINIMIZE  <<<<<<<<<< 

OBJ..  COST  =E=  SUM(I,  RD(I))  +  SUM(I,  BL(I))  ; 


>>>>>>>>>>  SUBJECT  TO  <<<<<<<<<< 


DREDl(I).. 

RD(I)  =G  = 

DRED2(I).. 

RD(I)  =G  = 

DBLUEl(I).. 

BL(I)  =G  = 

DBLUE2(I).. 

BL(I)  =G  = 

POWER((CRED(I)  -  ((A*BLUE(I)  +  B))),2); 
POWER((CRED(I)  -  ((A*BLUE(I)  +  B))),2); 
POWER((CBLUE(I)  -  ((C*RED(I)  +D))),2); 
POWER((CBLUE(I)  -  ((C*RED(I)  +D))),2); 
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MODEL  LSQ1  /OBJ,DRED  1,DBLUE  1  / ; 

SOLVE  LSQ1  USING  NLP  MINIMIZING  COST; 

A.FX  =  A.L; 

C.FX  =  C.L; 

MODEL  LSQ2  /OBJ,DRED2,DBLUE2/; 

SOLVE  LSQ2  USING  NLP  MINIMIZING  COST; 

- REPORTS - 

DISPLAY  COST.L,A  L,B.L,C.L,D.L 
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APPENDIX  E 


PROBABILITY  SURFACE  COMPARISON  CODE 

PROGRAM  datmat 

INTEGER  first(2),next(2. 100) 

REAL  ecdf(2,100) 

COMMON /list /first, next, ecdf 


INTEGER  rl(20),r2(20),lower,higher 

REAL  dmatrx(6, 6,500), cdfdat(6, 20), r(40) 

REAL  plow,phigh,tdiff(6,6) 

CHARACTER  *12  cfile 
DOUBLE  PRECISION  dseed 


C  dseed  =  7561883.d0 
dseed  =  829847.d0 

OPEN  (UNIT  =  8,FILE = ’cdfile.dat’, STATUS  =  ’old’) 
n  =  0 

111  CONTINUE 

n  =  n  +  1 

READ(8,21 1,END  =  55)  cfile 
OPEN(UNIT  =  9,FILE  =  cfile,STATUS  =  ’old’) 
m  =  0 

112  CONTINUE 

m  =  m  +  1 

READ(9,2 12,END  =  56)  cdfdat(n,m) 

GOTO  112 
56  CONTINUE 
CLOSE(9) 

GOTO  111 
55  CONTINUE 
CLOSE(8) 

OPEN(UNIT  =  10, FILE  =  ’datmat.out’,  STATUS  =  ’unknown’) 
OPEN(UNIT  =  1 1 , FILE  =  ’pdat.out’,STATUS  =  ’unknown’) 
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DO  113  i  =  1  ,  5 
DO  114  j  =  i+1  ,  6 
lower  =  0 
higher  =  0 

PRINT*,  ’Comparing  surfaces  for  :  \i,’  vs  \j 
DO  110  1  =  1  ,  20 
ecdf(l,l)  =  cdfdat(i,l) 
ecdf(2,l)  =  cdfdat(j,l) 

110  CONTINUE 

CALL  tcdf(tdiff(i,j)) 

WRITE(11,*)  ’i,j’,i,j,’  tdiff  =  ’,tdiff(i,j) 

PRINT  Vtdiff  =  ’,tdiff(i,j) 

DO  99  1  =  1  ,  20 
r(l)  =  cdfdat(i,l) 
r(l  +  20)  =  cdfdat(j,l) 

99  CONTINUE 

DO  116  1  =  1  ,  50 
CALL  roll(20,40,rl,r2,dseed) 

DO  115  k  =  1  ,  20 
ecdf(l,k)  =  r(rl(k)) 
ecdf(2,k)  =  r(r2(k)) 

1 15  CONTINUE 

CALL  tcdf(dmatrx(i,j,l)) 

IF(dmatrx(i,j,l)  .LT.  tdiff(i,j))  lower  =  lower  + 1 
IF(dmatrx(i,j,l)  .GT.  tdiff(i,j))  higher=  higher+1 

116  CONTINUE 

plow  =  FLOAT(lower)/50.0 
phigh  =  FLOAT(higher)/50.0 
114  CONTINUE 

113  CONTINUE 

DO  117  i  =  1  ,  5 
DO  118  j  =  i+1  ,  6 
DO  119  k  =  1  ,  50 
WRITE(  10,213)  dmatrx(i,j  k) 

1 19  CONTINUE 

118  CONTINUE 

117  CONTINUE 

211  FORMAT(A12) 

212  FORMAT(F10.7) 

213  FORMAT(F10.4) 

END 
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SUBROUTINE  tcdf(diff) 


INTEGER  first(2),next(2,100) 

REAL  ecdf(2,100) 

COMMON /list/first,next,ecdf 

REAL  deltax,deltay,x,y,diff, ratio 

DO  33  nc  =  1  ,2 
DO  32  n  =  1  ,  10 
next(nc,n)  =  0 
CALL  ordcdf(nc,n) 

32  CONTINUE 

33  CONTINUE 
diff  =  0.0 
deltax  =  1.0 
deltay  =  1.0 

DO  21  x  =  50,  200,  deltax 
Do  22  y  =  50,  200,  deltay 
ratio  =  (x/y)**2 
CALL  cdf(l, ratio, probl) 

CALL  cdf(2,ratio,prob2) 

diff  =  diff  +  (ab$(probl-prob2)  *  (dehax*  deltay)) 
22  CONTINUE 
21  CONTINUE 

100  FORMAT(F10.7) 

END 

SUBROUTINE  ordcdf(nc,n) 

INTEGER  count,test,ltest,ttest,maxcnt 

INTEGER  first(2),next(2,100) 

REAL  ecdf(2,100) 

COMMON /list/first,next,erdf 

maxcnt  =  100 
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IF  (n  .EQ.  1)  THEN 
first(nc)  =  1 
next(nc,l)  =  0 
count  =  1 
RETURN 

ELSEIF  (count  .EQ.  maxcnt)  THEN 
PRINT  *,  OVERFLOW  ERROR  >  >  >  CDF  LIST  IS  FULL  <  <  <’ 
STOP 
ELSE 

count  =  count  +  1 
test  =  first(nc) 

201  CONTINUE 

IF  (ecdf(nc,n)  .LT.  ecdf(nc.test))  THEN 
IF  (test  .EQ.  first(nc))  THEN 
first(nc)  =  n 
ELSE 

next(ncjtest)  =  n 
ENDIF 

next(nc,n)  =  test 
RETURN 
ELSE 

ttest  =  next(nc,test) 

IF  (ttest  .EQ.  0)  THEN 
next(nc,n)  =  0 
next(nc,test)  =  n 
RETURN 
ELSE 
ltest  =  test 
test  =  ttest 
GOTO  201 
ENDIF 
ENDIF 
ENDIF 
RETURN 

END 
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SUBROUTINE  cdf(nc,pt,prob) 

INTEGER  first(2),next(2,100) 
REAL  ecdf(2,100) 
COMMON/list/first,next,ecdf 

INTEGER  nc,test 
REAL  pt,prob 

num  =  0 
test  =  first(nc) 

400  CONTINUE 

IF  (pt  .LE.  ecdf(nc,test))  THEN 
prob  =  FLO  AT(  num)/ 10.0 
RETURN 
ELSE 

num  =  num  +  1 
test  =  next(nc,test) 

IF  (test  .EQ.  0)  THEN 
prob  =  1.0 
RETURN 
ENDIF 
ENDIF 
GOTO  400 
END 

SUBROUTINE  roll(n,m,r,rc,dseed) 

INTEGER  r(n),rc(n) 

DOUBLE  PRECISION  dseed 

DO  2  i  =  1  ,  n 
r(i)  =  0 

3  CONTINUE 

CALL  unif(dseed,u) 
num  =  INT((u*m)  +  .5) 

DO  1  j  =  1  ,  i 
IF  (r(j)  .EQ.  num)  GOTO  3 

1  CONTINUE 
r(i)  =  num 

2  CONTINUE 
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k  =  1 

DO  6  i  =  1  ,  m 
DO  4  j  =  1  ,  n 
IF  (r(j)  .EQ.  i)  GOTO  5 

4  CONTINUE 
rc(k)  =  i 

k  =  k  +  1 

5  CONTINUE 

6  CONTINUE 
RETURN 
END 

SUBROUTINE  UNIF(DSEED,U) 

C 

INTEGER  I 
REAL  U 

DOUBLE  PRECISION  DENOM.DSEED 
C 

DATA  DENOM /2 147483647.D0/ 

C 

DSEED  =  DMOD(  16807.D0*DSEED, DENOM) 
U  =  DSEED/DENOM 
C 

RETURN 

END 
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